'''
Created on Nov 5, 2009

@author: mkiyer
'''

class GTF2Feature(object):
    '''
    GTF2 file format support
    '''
    undefined_str = '.'
    attr_sep = ';'
    attr_tagvalue_sep = '='    
    attr_tags = ['ID', 'Name', 'Alias', 'Parent', 'Target', 'Gap', 
                 'Dbxref', 'Ontology_term']

    def __init__(self, args=None):
        self.seqid = self.undefined_str
        self.source = self.undefined_str
        self.type = self.undefined_str
        self.start = None
        self.end = None
        self.score = None
        self.strand = self.undefined_str
        self.phase = None
        self.attributes = self.undefined_str        
        if args is not None and isinstance(args, str):
            self._parse(args)
    
    def __str__(self):
        common_str = '\t'.join([self.seqid,
                                self.source,
                                self.type,
                                '.' if self.start is None else str(self.start),
                                '.' if self.end is None else str(self.end),
                                '.' if self.score is None else str(self.score),
                                self.strand,
                                '.' if self.phase is None else str(self.phase)])
        attr_fields = []
        for tag in self.attr_tags:
            if tag in self.attributes:
                attr_fields.append('%s%s%s' % (tag, self.attr_tagvalue_sep, self.attributes[tag]))
        return '\t'.join([common_str, self.attr_sep.join(attr_fields)])

    def _parse(self, input_string):
        '''
        Parse a string into a gff3 feature
        
        TODO: error checking during parsing, raise exceptions, etc 
        '''
        thisfields = input_string.split('\t')
        self.seqid = thisfields[0]
        self.source = thisfields[1]
        self.type = thisfields[2]
        self.start = int(thisfields[3]) if thisfields[3] != self.undefined_str else None
        self.end = int(thisfields[4]) if thisfields[4] != self.undefined_str else None            
        self.score = float(thisfields[5]) if thisfields[5] != self.undefined_str else None            
        self.strand = thisfields[6]
        self.phase = int(thisfields[7]) if thisfields[7] != self.undefined_str else None        
        if thisfields[8] != self.undefined_str:
            attrs = {}
            attr_strings = thisfields[8].split(self.attr_sep)
            for a in attr_strings:
                tag, value = a.split(self.attr_tagvalue_sep)
                attrs[tag] = value
            self.attributes = attrs

def read_gtf2_file(fhd, attr_defs=None):
    '''
    1. seqname - The name of the sequence. Must be a chromosome or scaffold.
    2. source - The program that generated this feature.
    3. feature - The name of this type of feature. Some examples of standard feature types are "CDS", "start_codon", "stop_codon", and "exon".
    4. start - The starting position of the feature in the sequence. The first base is numbered 1.
    5. end - The ending position of the feature (inclusive).
    6. score - A score between 0 and 1000. If the track line useScore attribute is set to 1 for this annotation data set, the score value will determine the level of gray in which this feature is displayed (higher numbers = darker gray). If there is no score value, enter ".".
    7. strand - Valid entries include '+', '-', or '.' (for don't know/don't care).
    8. frame - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'.

    chr1    Cufflinks       transcript      3109    3184    1000    .       .       gene_id "3064YAAXX_2.1"; transcript_id "3064YAAXX_2.1.0"; RPKM "0.1802096261"; frac "1.000000"; conf_lo "0.000000"; conf_hi "0.000000"; cov "0.108248";
    '''
    undefined_str = '.'
    attr_sep = ';'
    attr_tagvalue_sep = ' '    
    if attr_defs is None:
        attr_defs = {}
        
    for line in fhd:
        # read the GTF line
        line = line.strip()
        if line is None:
            continue
        if line.startswith('#'):
            continue
        fields = line.split('\t')
        seqname = fields[0]
        source = fields[1]
        feature = fields[2]
        start = int(fields[3])
        end = int(fields[4])
        score = float(fields[5])
        strand = fields[6]
        if not (strand == '+' or strand == '-'):
            strand = None
        frame = fields[7]
        if fields[8] != undefined_str:
            attrs = {}
            attr_strings = fields[8].split(attr_sep)
            for a in attr_strings:
                a = a.strip()
                if len(a) == 0:
                    continue
                tag, value = a.split(attr_tagvalue_sep)
                # remove quotes
                value = value.split('"')[1]
                # apply parsing function
                value = attr_defs[tag](value) if (tag in attr_defs and attr_defs[tag] is not None) else value
                attrs[tag] = value
        # convert to GFF3 style
        f = type('GFF3', (object,), dict())
        f.seqid = seqname
        f.source = source
        f.type = feature
        f.start = start
        f.end = end
        f.score = score
        f.strand = strand
        f.phase = frame
        f.attrs = attrs
        yield f


class Feature(object):
    def __init__(self, id):
        self.id = id
        self.type = self._feature_type
        self.parents = []
        self.children = []

    def add_parent(self, f):
        self.parents.append(f)
        f.children.append(self)

    def add_child(self, f):
        self.children.append(f)
        f.parents.append(self)
        
class Gene(Feature):
    _feature_type = 'gene'

    def __init__(self, id):
        super(Gene, self).__init__(id)
    
    @property
    def transcripts(self):
        return list(self.children)
            
class Transcript(Feature):
    _feature_type = 'transcript'
    
    def __init__(self, id, gene):
        super(Transcript, self).__init__(id)
        self.add_parent(gene)

    @property
    def gene(self):
        return list(self.parents)
    @property
    def exons(self):
        return list(self.children)    
    @property
    def chrom(self):
        return list(set([e.chrom for e in self.exons]))    
    @property
    def start(self):
        return min([e.start for e in self.exons])
    @property
    def end(self):
        return max([e.end for e in self.exons])
    @property
    def strand(self):
        return list(set([e.strand for e in self.exons]))
        

class Exon(Feature):
    _feature_type = 'exon'
    
    id_prefix = 'exon'
    id_count = 0
    
    @staticmethod
    def __get_id():
        id = Exon.id_prefix + str(Exon.id_count)
        Exon.id_count += 1
        return id
    
    def __init__(self, transcript, chrom, start, end, strand):
        id = self.__get_id()
        super(Exon, self).__init__(id)
        self.chrom = chrom
        self.start = start
        self.end = end
        self.strand = strand
        self.add_parent(transcript)
    
    @property
    def transcripts(self):
        return list(self.parents)

    def merge(self, other):
        merged_exon = copy.copy(self)
        # merge the exon intervals to the largest interval
        merged_exon.start = min(self.start, other.start)
        merged_exon.end = max(self.end, other.end)
        return merged_exon
        
class GeneExpr(Gene):
    id_prefix = 'gene'
    id_count = 0

    @staticmethod
    def __get_id():
        id = GeneExpr.id_prefix + str(GeneExpr.id_count)
        GeneExpr.id_count += 1
        return id
    
    def __init__(self, library_id):
        id = self.__get_id()
        super(GeneExpr, self).__init__(id)
        self.library_ids = [library_id]

    def add_library(self, library_id):
        self.library_ids.append(library_id)

class TranscriptExpr(Transcript):
    id_prefix = 'transcript'
    id_count = 0

    @staticmethod
    def __get_id():
        id = TranscriptExpr.id_prefix + str(TranscriptExpr.id_count)
        TranscriptExpr.id_count += 1
        return id
    
    def __init__(self, gene, library_id, rpkm):
        id = self.__get_id()
        super(TranscriptExpr, self).__init__(id, gene)
        self.library_ids = [library_id]
        self.covs = [rpkm]

    def add_library(self, library_id, rpkm):
        self.library_ids.append(library_id)
        self.covs.append(rpkm)

class ExonExpr(Exon):

    def __init__(self, transcript, chrom, start, end, strand, library_id, rpkm):
        super(ExonExpr, self).__init__(transcript, chrom, start, end, strand)
        self.library_ids = [library_id]
        self.covs = [rpkm]

    def add_library(self, library_id, rpkm):
        self.library_ids.append(library_id)
        self.covs.append(rpkm)
        
    def merge(self, other):
        merged_exon = copy.copy(self)
        # merge the exon intervals to the largest interval
        merged_exon.start = min(self.start, other.start)
        merged_exon.end = max(self.end, other.end)
        merged_exon.library_ids.extend(other.library_ids)
        merged_exon.covs.extend(other.covs)
        return merged_exon


cufflinks_attrs = {'gene_id': str,
                   'transcript_id': str,
                   'RPKM': float,
                   'frac': float,
                   'conf_lo': float,
                   'conf_hi': float,
                   'cov': float}

def build_gene_model(fhd):
    '''
    1. seqname - The name of the sequence. Must be a chromosome or scaffold.
    2. source - The program that generated this feature.
    3. feature - The name of this type of feature. Some examples of standard feature types are "CDS", "start_codon", "stop_codon", and "exon".
    4. start - The starting position of the feature in the sequence. The first base is numbered 1.
    5. end - The ending position of the feature (inclusive).
    6. score - A score between 0 and 1000. If the track line useScore attribute is set to 1 for this annotation data set, the score value will determine the level of gray in which this feature is displayed (higher numbers = darker gray). If there is no score value, enter ".".
    7. strand - Valid entries include '+', '-', or '.' (for don't know/don't care).
    8. frame - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'.

    chr1    Cufflinks       transcript      3109    3184    1000    .       .       gene_id "3064YAAXX_2.1"; transcript_id "3064YAAXX_2.1.0"; RPKM "0.1802096261"; frac "1.000000"; conf_lo "0.000000"; conf_hi "0.000000"; cov "0.108248";
    '''
    genes = {}
    transcripts = {}
    exons = {}
    
    linenum = 0
    debug_every = 10000
    
    for linenum, feature in enumerate(read_gtf2_file(fhd, attr_defs=cufflinks_attrs)):
        if linenum % debug_every == 0:
            print linenum
        # unpack feature
        chrom, start, end = feature.seqid, feature.start, feature.end
        type = feature.type
        strand = feature.strand
        attrs = feature.attrs

        if chrom != 'chr1':
            continue

        # update the gene model from this GTF line
        gene_id = attrs['gene_id']
        transcript_id = attrs['transcript_id']
        rpkm = attrs['RPKM']
        # create the gene if it does not already exist
        if gene_id in genes:
            g = genes[gene_id]
        else:
            # create new gene
            g = GeneExpr(gene_id)
            genes[gene_id] = g

        if type == 'transcript':
            # create new transcript of this gene
            t = TranscriptExpr(g, transcript_id, rpkm)                
            transcripts[transcript_id] = t            
        elif type == 'exon':
            if transcript_id in transcripts:
                t = transcripts[transcript_id]
            else:
                t = TranscriptExpr(g, transcript_id, rpkm)
                transcripts[transcript_id] = t            
            e = ExonExpr(t, chrom, start, end, strand, transcript_id, rpkm)
            exons[e.id] = e
    return genes, transcripts, exons

from bx.intervals.intersection import Interval, IntervalTree
import collections
import copy

def create_exome_tree(exons):
    trees = collections.defaultdict(lambda: IntervalTree())
    for e in exons:
        if e.chrom != 'chr1':
            continue
        # add a visited flag for exon merging
        e.visited = False
        trees[e.chrom].insert_interval(Interval(e.start, e.end, value=e))
    return trees

def get_overlap_percent(x1, x2, y1, y2):
    ostart, oend = sorted((x1, x2, y1, y2))[1:3]
    overlap_percent = min(float(oend - ostart)/(x2 - x1),
                          float(oend - ostart)/(y2 - y1))
    return overlap_percent

def compare_strands(strand1, strand2):
    if strand1 == '+' and strand2 == '-':
        return False
    elif strand1 == '-' and strand2 == '+':
        return False
    return True

def merge_exomes(e1, e2):    
    print 'building interval tree'
    trees = create_exome_tree(e1.itervalues())
    print 'merging'

    overlap_threshold = 0.90
    new_exome = {}
    multiple_hits = collections.defaultdict(lambda: 0)
    merged_exons = 0
    
    for exon in e2.itervalues():
        # search for intersection
        hits = trees[exon.chrom].find(exon.start, exon.end)
        if len(hits) == 0:
            # this exon has zero overlap, so just copy it into the new exome
            new_exome[exon.id] = exon
            continue

        if len(hits) > 1:
            # keep stats on multiple hits
            multiple_hits[len(hits)] += 1
            print 'multiple', exon.chrom, exon.start, exon.end, 'vs'
            for hit in hits:
                print '\t', 'multiple', hit.value.chrom, hit.value.start, hit.value.end 

        for hit in hits:
            matched_exon = hit.value
            # mark the matched exon as 'visited'
            matched_exon.visited = True
            # check that exons are on the same strand
            if not compare_strands(exon.strand, matched_exon.strand):
                print 'mismatched strands', exon.chrom, exon.start, exon.end, matched_exon.chrom, matched_exon.start, matched_exon.end
                # just copy the exons directly
                new_exome[exon.id] = exon
                new_exome[matched_exon.id] = matched_exon
                continue
            # calculate the % overlap between exons
            overlap = get_overlap_percent(exon.start, exon.end,
                                          matched_exon.start, matched_exon.end)
            # exon did not overlap enough, so they are
            # probably separate exons
            if overlap < overlap_threshold:
                print 'poor overlap', overlap, exon.chrom, exon.start, exon.end, matched_exon.chrom, matched_exon.start, matched_exon.end       
                # just copy the exon directly
                new_exome[exon.id] = exon
                new_exome[matched_exon.id] = matched_exon                
                continue
            # merge the two exons
            merged_exon = exon.merge(matched_exon)
            new_exome[merged_exon.id] = merged_exon
            merged_exons += 1
            print 'merged:', exon.chrom, exon.start, exon.end, '|', matched_exon.chrom, matched_exon.start, matched_exon.end
            print 'new exon:', merged_exon.chrom, merged_exon.start, merged_exon.end
            print 'new exon:', zip(merged_exon.library_ids, merged_exon.covs)

    print multiple_hits
    return new_exome

def merge_gene_models(g1, t1, e1, g2, t2, e2):
    print 'building interval tree'
    trees = create_exome_tree(e1.itervalues())
    print 'merging'
    
    new_genes, new_transcripts, new_exons = {}, {}, {}
    overlap_threshold = 0.9
    multiple_hits = 0

    for gene in g2.itervalues():
        
        # copy the gene and add to the new model
        new_g = copy.copy(gene)
        new_genes[new_g.id] = new_g
        
        for tx in gene.get_transcripts():            
            # prepare to match all the transcript's exons
            tx_matched = True
            matched_transcripts = None
            # preemptively build merged exons in case there is a match
            merged_exons = []
            # search exon by exon for overlapping exons between gene models
            for exon in tx.exons:
                hits = trees[exon.chrom].find(exon.start, exon.end)
                # this exon has zero overlap, so the transcript does not match
                if len(hits) == 0:
                    tx_matched = False
                    break
                # the exon overlaps greater than one reference exon... yikes
                if len(hits) > 1:
                    # TODO: multiple hits
                    multiple_hits += 1
                    for hit in hits:
                        print exon.chrom, exon.start, exon.end, 'vs', hit.value.chrom, hit.value.start, hit.value.end 
                    tx_matched = False
                    break                

                # calculate the % overlap between exons
                matched_exon = hits[0].value
                overlap = get_overlap_percent(exon.start, exon.end,
                                              matched_exon.start, matched_exon.end)
                if overlap < overlap_threshold:
                    tx_matched = False
                    break
                # at this point we know that this exon has matched
                # find all the transcripts that matched
                if matched_transcripts is None:
                    matched_transcripts = set(matched_exon.transcripts)
                else:
                    matched_transcripts.intersection_update(matched_exon.transcripts)
                # add the union of both exons to the potential new transcript
                merged_exons.append((min(exon.start, matched_exon.start),
                                     max(exon.end, matched_exon.end)))

            if tx_matched == False:
                # create the new transcript
                new_t = copy.copy(tx)
            else:
                # deep copy the transcript
                new_t = copy.deepcopy(tx)
                
            new_transcripts[new_t.id] = new_t
            
            # merge if the transcript matched
            if tx_matched == True:
                
                new_g = copy.copy(g2)
                
                t2_gene_id = tx.get_parent().id
                # create the gene if it does not already exist
                if t2_gene in new_genes:
                    new_g = genes[t2_gene_id]
                else:
                    # create new gene
                    new_g = GeneExpr(t2_gene_id)
                    new_genes[t2_gene_id] = new_g
    
            if type == 'transcript':
                # create new transcript of this gene
                t = TranscriptExpr(g, chrom, start, end, strand, transcript_id, rpkm)                
                transcripts[transcript_id] = t            
            elif type == 'exon':
                if transcript_id in transcripts:
                    t = transcripts[transcript_id]
                else:
                    t = TranscriptExpr(g, chrom, start, end, strand, transcript_id, rpkm)
                    transcripts[transcript_id] = t            
                e = ExonExpr(t, chrom, start, end, transcript_id, rpkm)
                exons[e.id] = e


    matched_coverage = []    
    # intersect each interval
    for a1 in annots:
        id1 = a1.name.split('.')[0] 
        hits = trees[a1.chrom].find(a1.start, a1.end)
        true_hits = []
        library_best_hits = {}
        library_best_overlap = collections.defaultdict(lambda: 0.0)
                
        for hit in hits:
            a2 = hit.value
            id2 = a2.name.split('.')[0]                
            # ignore hits from same library (isoforms)
            if id1 == id2:
                continue
            # find overlapping interval
            ostart, oend = sorted((a1.start, a1.end, a2.start, a2.end))[1:3]
            overlap_percent = max(float(oend - ostart)/(a1.end - a1.start),
                                  float(oend - ostart)/(a2.end - a2.start))
            # TODO: set an overlap threshold?
            #library_hits[id2].append((a2, overlap_percent))
            # choose the isoform with best overlap for each library
            if overlap_percent > library_best_overlap[id2]:
                library_best_overlap[id2] = overlap_percent
                library_best_hits[id2] = a2
        # organize hits by sample
        sample_coverages = collections.defaultdict(lambda: [])
        for transcript in library_best_hits.values():
            sample_coverages[transcript.sample].append(transcript.score)
        # TODO: how to find average?  with or without "zeros"?
        # for now, do this without the zeros included, assuming the lanes
        # are poor quality and depth of coverage is poor
        avg1 = np.average(sample_coverages[sample1])
        avg2 = np.average(sample_coverages[sample2])
        matched_coverage.append((avg1, avg2))        
        if avg1 == 0: avg1 = 0.01
        if avg2 == 0: avg2 = 0.01
        if max(avg1, avg2) > 10.0 and abs(math.log(avg1/avg2)) >= 1.0:
            print '\t'.join([a1.chrom, str(a1.start), str(a1.end), str(avg1), str(avg2)])
        #    print a1.chrom, a1.start, a1.end, a1.name, a1.score, a2.chrom, a2.start, a2.end, a2.name, a2.score        
        

def parse_transcripts_tmap(fhd):
    '''
    ref_gene_id     ref_id  class_code      cuff_gene_id    cuff_id FMI     RPKM    <unknown>    <unknown>    cov     len     major_iso_id
    ['OTTHUMT00000002844', 'OTTHUMT00000002844', 'c', '3064YAAXX_2.1', '3064YAAXX_2.1.0', '100', '0.164120', '0.000000', '0.000000', '0.108248', '76', '3064YAAXX_2.1.0']
    '''
    # first line is header
    fhd.readline()
    for thisline in fhd:        
        if not thisline:
            continue
        if thisline.startswith('#'):
            continue
        thisline = thisline.strip()
        if not thisline:
            continue
        thisfields = thisline.split('\t')
        t = type('TranscriptMap', (object,), dict())
        t.line = thisline
        t.ref_gene_id = thisfields[0] if thisfields[0] is not '-' else None
        t.ref_id = thisfields[1] if thisfields[1] is not '-' else None
        t.class_code = thisfields[2]
        t.cuff_gene_id = thisfields[3]
        t.cuff_id = thisfields[4]
        t.fmi = int(thisfields[5])
        t.rpkm = float(thisfields[6])
        t.cov = float(thisfields[9])
        t.length = int(thisfields[10])
        t.major_iso_id = thisfields[11]
        yield t

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import math

def scatter_hist(x, y):
    fig = plt.figure(1, figsize=(5.5,5.5))    
    from mpl_toolkits.axes_grid import make_axes_locatable
    
    axScatter = plt.subplot(111)
    divider = make_axes_locatable(axScatter)
    
    # create a new axes with a height of 1.2 inch above the axScatter
    axHistx = divider.new_vertical(1.2, pad=0.1, sharex=axScatter)
    
    # create a new axes with a width of 1.2 inch on the right side of the
    # axScatter
    axHisty = divider.new_horizontal(1.2, pad=0.1, sharey=axScatter)
    
    fig.add_axes(axHistx)
    fig.add_axes(axHisty)
    
    # make some labels invisible
    plt.setp(axHistx.get_xticklabels() + axHisty.get_yticklabels(),
             visible=False)
    
    # the scatter plot:
    axScatter.scatter(x, y)
    axScatter.set_aspect(1.)
    
    # now determine nice limits by hand:
    #binwidth = 0.25
    binwidth = 0.01
    xymax = np.max( [np.max(np.fabs(x)), np.max(np.fabs(y))] )
    lim = ( int(xymax/binwidth) + 1) * binwidth
    
    bins = np.arange(0, lim + binwidth, binwidth)
    axHistx.hist(x, bins=bins)
    axHisty.hist(y, bins=bins, orientation='horizontal')
    
    # the xaxis of axHistx and yaxis of axHisty are shared with axScatter,
    # thus there is no need to manually adjust the xlim and ylim of these
    # axis.
    
    #axHistx.axis["bottom"].major_ticklabels.set_visible(False)
    for tl in axHistx.get_xticklabels():
        tl.set_visible(False)
    axHistx.set_yticks([0, 50, 100])
    
    #axHisty.axis["left"].major_ticklabels.set_visible(False)
    for tl in axHisty.get_yticklabels():
        tl.set_visible(False)
    axHisty.set_xticks([0, 50, 100])
    
    plt.draw()
    plt.show()
    #plt.savefig("a.pdf")
        
def plot_scatter(x, y, 
                 outfile=None,
                 xlabel='Control Group',
                 ylabel='Experimental Group',
                 title='Coverage scatter plot'):  
    fig = plt.figure()
    ax = fig.add_subplot(111)
    plt.title(title)

    #sorted_data = sorted(data, operator.itemgetter(0))
    xdata = x
    ydata = y
    #xdata = [d[0] for d in data]
    #ydata = [d[1] for d in data]

    # size?
    # s = 20*np.log(np.array([x.ratio for x in intervals]))        
    # color?
    # c = [x.category * cm.jet.N / len(category_names)] * len(intervals)
    # c = x.category * cm.jet.N / len(category_names)
    plt.scatter(xdata, ydata, 
                s=30, c='b', marker='o', cmap=None, norm=None,
                vmin=None, vmax=None, alpha=1.0, linewidths=None,
                verts=None)

    if outfile is None:
        plt.show()
    else:
        plt.savefig(outfile)
    plt.close()
    
    return 

import sys
import array

if __name__ == '__main__':
    
    gtf_file = open(sys.argv[1])    
    print 'reading first file'
    g1, t1, exome = build_gene_model(gtf_file)
    gtf_file.close()
    print len(g1), 'genes'
    print len(t1), 'transcripts'
    print len(exome), 'exons'

    for gtf_file in sys.argv[2:]:
        print 'reading file'
        g2, t2, e2 = build_gene_model(open(gtf_file))
        print len(g2), 'genes'
        print len(t2), 'transcripts'
        print len(e2), 'exons'
        
        print 'merging'
        exome = merge_exomes(exome, e2)
        print len(exome), 'exons'
        print len(exome), 'exons'
        print len(exome), 'merged exome'

    sys.exit(0)
    
    '''
    1      =      Match
    2     c     Contained     
    3     j     New isoform     
    4     e     A single exon transcript overlapping a reference exon and at least 10 bp of a reference intron, indicating a possible pre-mRNA fragment.     
    5     i     A single exon transcript falling entirely with a reference intron     
    6     r     Repeat. Currently determined by looking at the reference sequence and applied to transcripts where at least 50% of the bases are lower case     
    7     u     Unknown, intergenic transcript
    '''
    exprs = array.array('f')
    lengths = array.array('f')
    
    novel_class_codes = set(['i', 'u', 'j'])
     
    for tmap in parse_transcripts_tmap(open(sys.argv[1])):
        if tmap.length < 500:
            if tmap.rpkm == 0.0:
                tmap.rpkm = 0.001 
            exprs.append(math.log(tmap.rpkm))
            lengths.append(tmap.length)
    plot_scatter(lengths, exprs)
    #scatter_hist(lengths, exprs)
